﻿#include "MultiSystem.h"
#include "ForwardStep.h"
#include "IsoVortex.h"
#include "OverlapMesh.h"
#include "PassFlow.h"
#include <vtkm/VectorAnalysis.h>

struct MeshModifier1 : vtkm::worklet::WorkletCellNeighborhood
{
  using ControlSignature = void(CellSetIn,
                                FieldInNeighborhood centroid,
                                FieldOut blk_id,
                                FieldOut bnd_id);
  using ExecutionSignature = void(Boundary, _2, _3, _4);

  template<typename NeighIn, typename BndIdType>
  VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary,
                            const NeighIn& centroid,
                            vtkm::Int16& blk_id,
                            BndIdType& bnd_id) const
  {
    blk_id = 0;
    if (!IsCellValid(centroid.Get(0, 0, 0)))
    {
      blk_id = -1; // 不可用单元
    }

    bnd_id = { -1, -1, -1, -1 };
    if (boundary.IJK[0] == 0) // left, inlet
    {
      bnd_id[3] = 3;
    }
    if (boundary.IJK[0] == boundary.PointDimensions[0] - 1) // right, outlet
    {
      bnd_id[1] = 1;
    }
    if (boundary.IJK[1] == boundary.PointDimensions[1] - 1) // top, wall
    {
      bnd_id[2] = 2;
    }
    if (boundary.IJK[1] == 0) // bottom
    {
      bnd_id[0] = 0;
    }

    if (!IsCellValid(centroid.Get(-1, 0, 0)))
    {
      bnd_id[3] = 4;
    }
    if (!IsCellValid(centroid.Get(1, 0, 0)))
    {
      bnd_id[1] = 4;
    }
    if (!IsCellValid(centroid.Get(0, -1, 0)))
    {
      bnd_id[0] = 4;
    }
    if (!IsCellValid(centroid.Get(0, 1, 0)))
    {
      bnd_id[2] = 4;
    }
  }

  VTKM_EXEC bool IsCellValid(Vec3f p) const
  {
    //return Square(p);
    return Sphere(p);
    //return Naca0012(p);
    //return ForwardStep(p);
  }

  VTKM_EXEC bool Square(const Vec3f& p) const
  {
    Real x = p[0];
    Real y = p[1];
    Real a = 0.1;
    if (x < a && x > -a && y < a && y > -a)
    {
      return false;
    }
    else if (x < 0.4 && x > 0.3 && y < 0.7 && y > 0.6)
    {
      return false;
    }
    else if (x < 0.8 && x > 0.6 && y < 0.4 && y > 0.2)
    {
      return false;
    }
    else if (x < 1.5 && x > 1.3 && y < 0.2 && y > -0.1)
    {
      return false;
    }
    else
    {
      return true;
    }
  }

  VTKM_EXEC bool Sphere(const Vec3f& p) const
  {
    if (vtkm::Magnitude(p) < 0.1)
    {
      return false;
    }
    /*if (vtkm::Magnitude(p - Vec3f{ 0.5, 0.2, 0 }) < 0.2)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.5, -0.2, 0 }) < 0.2)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.7, 0.2, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.7, -0.2, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.9, 0.3, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.9, -0.3, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 1.7, -0.3, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 1.5, 0.5, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 2.5, 0.3, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 2.5, -0.3, 0 }) < 0.3)
    {
      return false;
    }*/

    return true;
  }

  VTKM_EXEC bool Naca0012(const Vec3f& p) const
  {
    Real x = p[0];
    Real y = p[1];

    Real a0 = +0.2969, a1 = -0.1260, a2 = -0.3516, a3 = +0.2843, a4 = -0.1036;
    Real t = 0.12;

    if (x < 0 || x > 1)
    {
      return true;
    }
    else
    {
      Real yt =
        5 * t * (a0 * vtkm::Sqrt(x) + a1 * x + a2 * x * x + a3 * x * x * x + a4 * x * x * x * x);
      if (y > yt || y < -yt)
      {
        return true;
      }
      else
      {
        return false;
      }
    }
  }

  VTKM_EXEC bool ForwardStep(const Vec3f& p) const
  {
    vtkm::Bounds box(Vec3f{ 0.6, 0, 0 }, Vec3f{ 3, 0.2, 0 });
    if (box.Contains(p))
    {
      return false;
    }
    return true;
  }
};

struct MeshModifier2 : vtkm::worklet::WorkletCellNeighborhood
{
  using ControlSignature = void(CellSetIn,
                                FieldInNeighborhood centroid,
                                FieldOut blk_id,
                                FieldOut bnd_id);
  using ExecutionSignature = void(Boundary, _2, _3, _4);

  template<typename NeighIn, typename BndIdType>
  VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary,
                            const NeighIn& centroid,
                            vtkm::Int16& blk_id,
                            BndIdType& bnd_id) const
  {
    blk_id = 0;
    if (!IsCellValid(centroid.Get(0, 0, 0)))
    {
      blk_id = -1; // 不可用单元
    }

    bnd_id = { -1, -1, -1, -1 };
    if (boundary.IJK[0] == 0) // left, inlet
    {
      bnd_id[3] = 3;
    }
    if (boundary.IJK[0] == boundary.PointDimensions[0] - 1) // right, outlet
    {
      bnd_id[1] = 1;
    }
    if (boundary.IJK[1] == boundary.PointDimensions[1] - 1) // top, wall
    {
      bnd_id[2] = 2;
    }
    if (boundary.IJK[1] == 0) // bottom
    {
      bnd_id[0] = 0;
    }

    if (!IsCellValid(centroid.Get(-1, 0, 0)))
    {
      bnd_id[3] = 11;
    }
    if (!IsCellValid(centroid.Get(1, 0, 0)))
    {
      bnd_id[1] = 11;
    }
    if (!IsCellValid(centroid.Get(0, -1, 0)))
    {
      bnd_id[0] = 11;
    }
    if (!IsCellValid(centroid.Get(0, 1, 0)))
    {
      bnd_id[2] = 11;
    }
  }

  VTKM_EXEC bool IsCellValid(Vec3f p) const
  {
    //return Square(p);
    return Sphere(p);
    //return Naca0012(p);
    //return ForwardStep(p);
  }

  VTKM_EXEC bool Square(const Vec3f& p) const
  {
    Real x = p[0];
    Real y = p[1];
    Real a = 0.1;
    if (x < a && x > -a && y < a && y > -a)
    {
      return false;
    }
    else if (x < 0.4 && x > 0.3 && y < 0.7 && y > 0.6)
    {
      return false;
    }
    else if (x < 0.8 && x > 0.6 && y < 0.4 && y > 0.2)
    {
      return false;
    }
    else if (x < 1.5 && x > 1.3 && y < 0.2 && y > -0.1)
    {
      return false;
    }
    else
    {
      return true;
    }
  }

  VTKM_EXEC bool Sphere(const Vec3f& p) const
  {
    if (vtkm::Magnitude(p) < 0.15)
    {
      return false;
    }
    return true;
  }

  VTKM_EXEC bool Naca0012(const Vec3f& p) const
  {
    Real x = p[0];
    Real y = p[1];

    Real a0 = +0.2969, a1 = -0.1260, a2 = -0.3516, a3 = +0.2843, a4 = -0.1036;
    Real t = 0.12;

    if (x < 0 || x > 1)
    {
      return true;
    }
    else
    {
      Real yt =
        5 * t * (a0 * vtkm::Sqrt(x) + a1 * x + a2 * x * x + a3 * x * x * x + a4 * x * x * x * x);
      if (y > yt || y < -yt)
      {
        return true;
      }
      else
      {
        return false;
      }
    }
  }

  VTKM_EXEC bool ForwardStep(const Vec3f& p) const
  {
    vtkm::Bounds box(Vec3f{ 0.6, 0, 0 }, Vec3f{ 3, 0.2, 0 });
    if (box.Contains(p))
    {
      return false;
    }
    return true;
  }
};

RegisterObject(MultiSystem);
MultiSystem::MultiSystem(const Configuration& cfg)
  : System(cfg)
{

  auto& overlap_mesh = static_cast<OverlapMesh&>(_mesh);

  {
    Configuration cfg1;
    cfg1.Add<Application*>("_app", &_app);
    Mesh* structured_mesh = overlap_mesh.GetStructuredMesh().get();
    cfg1.Add<Mesh*>("_mesh", structured_mesh);
    _structured_sys = std::make_shared<PassFlow>(cfg1);
  }

  {
    Configuration cfg1;
    cfg1.Add<Application*>("_app", &_app);
    Mesh* unstructured_mesh = overlap_mesh.GetUnStructuredMesh().get();
    cfg1.Add<Mesh*>("_mesh", unstructured_mesh);
    _unstructured_sys = std::make_shared<PassFlow>(cfg1);
  }

  _structured_sys->SetOther(_unstructured_sys);
  _unstructured_sys->SetOther(_structured_sys);
}

void MultiSystem::Init()
{
  _structured_sys->Init();
  _unstructured_sys->Init();

  {
    Mesh& mesh = _structured_sys->GetMesh();
    vtkm::cont::Invoker invoker;
    const auto& centroid = mesh.GetCellField("center");
    ArrayHandle<vtkm::Int16> blk_id;
    ArrayHandle<Vec<vtkm::Int16, 4>> bnd_id;

    vtkm::cont::Invoker{}(MeshModifier1{}, mesh.GetCellSet(), centroid, blk_id, bnd_id);

    mesh.GetDataSet().AddCellField("blk_id", blk_id);
    mesh.GetDataSet().AddCellField("bnd_id", bnd_id);
  }

  {
    Mesh& mesh = _unstructured_sys->GetMesh();
    vtkm::cont::Invoker invoker;
    const auto& centroid = mesh.GetCellField("center");
    ArrayHandle<vtkm::Int16> blk_id;
    ArrayHandle<Vec<vtkm::Int16, 4>> bnd_id;

    vtkm::cont::Invoker{}(MeshModifier2{}, mesh.GetCellSet(), centroid, blk_id, bnd_id);

    mesh.GetDataSet().AddCellField("blk_id", blk_id);
    mesh.GetDataSet().AddCellField("bnd_id", bnd_id);
  }

  console::Success("MultiSystem::Init 初始条件设置完成");
}

void MultiSystem::Evolve()
{
  _structured_sys->Evolve();
  _unstructured_sys->Evolve();


  console::Info("MultiSystem::Evolve运行完成，");
}

std::shared_ptr<PassFlow>& MultiSystem::GetStructuredSystem()
{
  return _structured_sys;
}

std::shared_ptr<PassFlow>& MultiSystem::GetUnStructuredSystem()
{
  return _unstructured_sys;
}
